
##### LOAD PACKAGES AND DATA #####
library(tidyverse)
library(haven)
library(MASS)
library(survey)
library(stargazer)
library(ggplot2)
library(lavaan)
library(interflex)

Lucid2020 <- read_stata("Lucid2020_Coded.dta")
Lucid2021 <- read_stata("Lucid2021_Coded.dta")
ANESGSS2020 <- read_stata("ANESGSS2020_Coded.dta")

#### post FUNCTION ####

setGeneric("post", 
           function(model,x1name=NULL,x1vals=NULL,x2name=NULL,x2vals=NULL,holds=NULL,seed=NULL,
                    n.sims=1000,cut=NULL,quantiles=c(.025,.975),did=NULL,weights=NULL, digits=2){
             standardGeneric("post")
           }
)

setClassUnion("arrayORNULL", c("array","NULL"))
setClassUnion("listORcharacter", c("list","character"))

setClass("post", 
         slots = c(est = "array", 
                   did = "arrayORNULL", 
                   sims = "array", 
                   model = "character", 
                   link = "listORcharacter", 
                   quantiles = "numeric", 
                   call = "call")
)


post.glm <- function(model,x1name=NULL,x1vals=NULL,x2name=NULL,x2vals=NULL,holds=NULL,seed=NULL,
                     n.sims=1000,cut=NULL,quantiles=c(.025,.975),did=NULL,weights=NULL, digits=2){
  
  call <- match.call()
  
  set.seed(seed)
  sims <- postSim(model, n.sims=n.sims)
  
  if (family(model)[2]=="identity"){link <- identity} 
  else if (family(model)[2]=="probit"){link <- pnorm}
  else if (family(model)[2]=="logit"){link <- plogis}
  else if (family(model)[2]=="log"){link <- exp}
  else if (family(model)[2]=="cloglog"){link <- function(x){1-exp(-exp(x))}}
  else {stop("Link function is not supported")}
  
  if (is.null(weights)){wi <- c(rep(1, length(model$model[,1])))} else{wi <- weights}
  n.obs <- length(model.matrix(model)[,1])
  k <- length(model.matrix(model)[1,])
  n.q <- length(quantiles)
  
  if (is.null(x1name)){
    X <- array(NA, c(n.obs,k))
    newdata <- data.frame(model$model)
    if (!is.null(holds)){
      for (i in 1:length(holds)){
        newdata[ ,names(holds)[i]] <- as.numeric(holds[i])
      }
    }
    X <- t(model.matrix(lm(formula(model), data=newdata)))
    l1 <- array(NA, c(nrow(sims@coef),1))
    l1[,1] <- apply(link(sims@coef %*% X), 1, function(x) weighted.mean(x, wi))
    l2 <- array(NA, c(1,n.q+1))
    l2[1,1] <- weighted.mean(link(coef(model) %*% X), wi)
    l2[1,2:(n.q+1)] <- quantile(l1, probs=quantiles)
    colnames(l2) <- c("est",quantiles)
    
    ans <- new("post", 
               est=round(l2, digits=digits),
               did=NULL,
               sims=l1, 
               model=class(model), 
               link=family(model)[2], 
               quantiles=quantiles, 
               call=call)
    return(ans)
  }
  
  else if (is.null(x2name)){
    
    n.x1 <- length(x1vals)
    X <- array(NA, c(n.obs,k,n.x1))
    
    for (i in 1:(n.x1)){
      newdata <- data.frame(model$model)
      if (!is.null(holds)){
        for (j in 1:length(holds)){
          newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
        }
      }
      newdata[ ,x1name] <- x1vals[i]
      X[ , ,i] <- model.matrix(lm(formula(model), data=newdata))
    }
    
    X <- aperm(X, c(2,1,3))
    l1 <- apply(apply(X, c(2,3), function(x) link(sims@coef %*% x)), c(1,3), function(x) weighted.mean(x, wi))
    l2 <- array(NA, c(n.x1+1,n.q+1))
    l2[1:n.x1,1] <- apply(apply(X, c(2,3), function(x) link(coef(model) %*% x)), 2, function(x) weighted.mean(x, wi))
    l2[1:n.x1,2:(n.q+1)] <- aperm(apply(l1, 2, function(x) quantile(x, probs=quantiles)))
    l2[nrow(l2),1] <- l2[n.x1,1] - l2[1,1]
    l2[nrow(l2),2:(n.q+1)] <- quantile(l1[ ,n.x1] - l1[ ,1], probs=quantiles)
    rownames(l2) <- c(paste(c(rep(paste(x1name,"="),n.x1),
                              paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),
                            c(x1vals,"")))  
    colnames(l2) <- c("est",quantiles)
    
    ans <- new("post", 
               est=round(l2, digits=digits), 
               did=NULL, 
               sims=l1, 
               model=class(model), 
               link=family(model)[2], 
               quantiles=quantiles, 
               call=call)
    return(ans)
  } 
  
  else{
    
    n.x1 <- length(x1vals)
    n.x2 <- length(x2vals)
    X <- array(NA, c(n.obs,k,n.x1,n.x2))
    
    for (j in 1:n.x2){
      for (i in 1:n.x1){
        newdata <- data.frame(model$model)
        if (!is.null(holds)){
          for (k in 1:length(holds)){
            newdata[ ,names(holds)[k]] <- as.numeric(holds[k])
          }
        }
        newdata[ ,x1name] <- x1vals[i]
        newdata[ ,x2name] <- x2vals[j]
        X[ , ,i,j] <- model.matrix(lm(formula(model), data=newdata))
      }
    }
    
    X <- aperm(X, c(2,1,3,4))
    l1 <- apply(apply(X, c(2,3,4), function(x) link(sims@coef %*% x)), c(1,3,4), function(x) weighted.mean(x, wi))
    l2 <- array(NA, c(n.x1+1,n.q+1,n.x2))
    l2[1:n.x1,1,1:n.x2] <- apply(apply(X, c(2,3,4), function(x) link(coef(model) %*% x)), c(2,3), function(x) weighted.mean(x, wi))
    l2[1:n.x1,2:(n.q+1),1:n.x2] <- aperm(apply(l1, c(2,3), function(x) quantile(x, probs=quantiles)), c(2,1,3))
    l2[nrow(l2),1,1:n.x2] <- l2[n.x1,1,1:n.x2] - l2[1,1,1:n.x2]
    l2[nrow(l2),2:(n.q+1),1:n.x2] <- apply(l1[ ,n.x1,1:n.x2] - l1[ ,1,1:n.x2], 2, function(x) quantile(x, probs=quantiles))
    dimnames(l2) <- list(paste(c(rep(paste(x1name,"="),n.x1),paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),c(x1vals,"")),
                         c("est",quantiles),
                         paste(c(rep(paste(x2name,"="),n.x2)),
                               c(x2vals)))
    
    if (is.null(did)){did <- c(x2vals[1],x2vals[n.x2])} else{did <- did}
    l3 <- array(NA, c(1,n.q+1))
    l3[1,1] <- (l2[n.x1,1,match(did[2],x2vals)] - l2[1,1,match(did[2],x2vals)]) - (l2[n.x1,1,match(did[1],x2vals)] - l2[1,1,match(did[1],x2vals)])
    l3[1,2:(n.q+1)] <- quantile((l1[ ,n.x1,match(did[2],x2vals)] - l1[ ,1,match(did[2],x2vals)]) -  (l1[ ,n.x1,match(did[1],x2vals)] - l1[ ,1,match(did[1],x2vals)]), probs=quantiles)
    dimnames(l3) <- list("did",c("est",quantiles)) 
    
    ans <- new("post", 
               est=round(l2, digits=digits), 
               did=round(l3, digits=digits), 
               sims=l1, 
               model=class(model), 
               link=family(model)[2], 
               quantiles=quantiles, 
               call=call)
    return(ans)
  }
}


post.polr <- function(model,x1name=NULL,x1vals=NULL,x2name=NULL,x2vals=NULL,holds=NULL,seed=NULL,
                      n.sims=1000,cut=NULL,quantiles=c(.025,.975),did=NULL,weights=NULL, digits=3){
  
  call <- match.call()
  
  set.seed(seed)
  sims <- suppressMessages(postSim(model, n.sims=n.sims))
  
  if (is.null(weights)){wi <- c(rep(1, length(model$model[,1])))} else{wi <- weights}
  
  if (model$method=="probit"){link <- pnorm}
  else if (model$method=="logistic"){link <- plogis}
  else if (model$method=="cloglog"){link <- function(x){1-exp(-exp(x))}}
  else {stop("Link function is not supported")}
  
  n.obs <- length(model$model[,1])
  k <- length(model.matrix(polr(getCall(model)$formula, model$model))[1,])
  n.q <- length(quantiles)
  n.y <- length(levels(model$model[,1]))
  n.z <- length(model$zeta)
  tau <- array(NA, c(n.sims,n.z+2))
  zeta <- array(NA, c(n.z+2))
  tau[,1] <- zeta[1] <- -Inf
  tau[,2:(ncol(tau)-1)] <- sims@zeta[,1:n.z]
  zeta[2:(ncol(tau)-1)] <- model$zeta[1:n.z]
  tau[,ncol(tau)] <- zeta[ncol(tau)] <- Inf
  beta <- sims@coef
  
  if (is.null(cut)){
    
    if (is.null(x1name)){
      
      X_temp <- array(NA, c(n.obs,k))
      X <- array(NA, c(n.obs,k-1))
      
      newdata <- data.frame(model$model)
      if (!is.null(holds)){
        for (j in 1:length(holds)){
          newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
        }
      }
      X_temp[ , ] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
      X[ , ] <- X_temp[,-1]
      X <- aperm(X)
      
      l1 <- array(NA, c(n.sims, n.obs, n.y))
      l2 <- array(NA, c(n.sims,n.y))
      l3 <- array(NA, c(n.y,n.q+1))
      for (z in 1:n.y){
        l1[,,z] <- link(tau[,z+1] - beta %*% X) - link(tau[,z] - beta %*% X)
        l2[,z] <- apply(l1[,,z], 1, function(x) weighted.mean(x, wi))
        l3[z,1] <- weighted.mean(link(zeta[z+1] - coef(model) %*% X) - link(zeta[z] - coef(model) %*% X), wi)
        l3[z,2:(n.q+1)] <- quantile(l2[,z], probs=quantiles)
      }
      
      rownames(l3) <- paste(c(rep("Y =",n.y)), c(1:n.y))
      colnames(l3) <- c("est",quantiles)
      
      ans <- new("post", 
                 est=round(l3, digits=digits), 
                 did=NULL, 
                 sims=l2, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
    
    else if (is.null(x2name)){
      
      n.x1 <- length(x1vals)
      
      X_temp <- array(NA, c(n.obs,k,n.x1))
      X <- array(NA, c(n.obs,k-1,n.x1))
      
      for (i in 1:(n.x1)){
        newdata <- data.frame(model$model)
        if (!is.null(holds)){
          for (j in 1:length(holds)){
            newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
          }
        }
        newdata[ ,x1name] <- x1vals[i]
        X_temp[ , ,i] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
        X[ , ,i] <- X_temp[,-1,i]
      }
      
      l1 <- array(NA, c(n.sims, n.obs, n.x1, n.y))
      l2 <- array(NA, c(n.sims, n.x1, n.y))
      l3 <- array(NA, c(n.x1+1, n.q+1, n.y))
      X <- aperm(X, c(2,1,3))
      for (z in 1:n.y){
        l1[,,,z] <- apply(X, c(2,3), function(x) (link(tau[,z+1] - beta %*% x) - link(tau[,z] - beta %*% x)))
        l2[,,z] <- apply(l1[,,,z], c(1,3), function(x) weighted.mean(x, wi))
        l3[1:n.x1,1,z] <- apply(apply(X, c(2,3), function(x) (link(zeta[z+1] - coef(model) %*% x) - link(zeta[z] - coef(model) %*% x))), 2, function(x) weighted.mean(x, wi))
        l3[1:n.x1,2:(n.q+1),z] <- t(apply(l2[,,z], 2, function(x) quantile(x, probs=quantiles)))
        l3[nrow(l3),1,z] <- l3[n.x1,1,z] - l3[1,1,z]
        l3[nrow(l3),2:(n.q+1),z] <- quantile(l2[ ,n.x1,z] - l2[ ,1,z], probs=quantiles)
      }
      
      dimnames(l3) <- list(paste(c(rep(paste(x1name,"="),n.x1),paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),c(x1vals,"")),
                           c("est",quantiles),
                           paste(c(rep("Y =",length(levels(model$model[,1])))),
                                 c(1:length(levels(model$model[,1]))))) 
      
      ans <- new("post", 
                 est=round(l3, digits=digits), 
                 did=NULL, 
                 sims=l2, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
    
    else{
      
      n.x1 <- length(x1vals)
      n.x2 <- length(x2vals)
      
      X_temp <- array(NA, c(n.obs,k,n.x1,n.x2))
      X <- array(NA, c(n.obs,k-1,n.x1,n.x2))
      
      for (j in 1:n.x2){
        for (i in 1:(n.x1)){
          newdata <- data.frame(model$model)
          if (!is.null(holds)){
            for (k in 1:length(holds)){
              newdata[ ,names(holds)[k]] <- as.numeric(holds[k])
            }
          }
          newdata[ ,x1name] <- x1vals[i]
          newdata[ ,x2name] <- x2vals[j]
          X_temp[ , ,i,j] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
          X[ , ,i,j] <- X_temp[,-1,i,j]
        }
      }
      
      X <- aperm(X, c(2,1,3,4))
      
      l1 <- array(NA, c(n.sims, n.obs, n.x1, n.x2, n.y))
      l2 <- array(NA, c(n.sims, n.x1, n.x2, n.y))
      l3 <- array(NA, c(n.x1+1, n.q+1, n.x2, n.y))
      for (z in 1:n.y){
        l1[,,,,z] <- apply(X, c(2,3,4), function(x) (link(tau[,z+1] - beta %*% x) - link(tau[,z] - beta %*% x)))
        l2[,,,z] <- apply(l1[,,,,z], c(1,3,4), function(x) weighted.mean(x, wi))
        l3[1:n.x1,1,1:n.x2,z] <- apply(apply(X, c(2,3,4), function(x) (link(zeta[z+1] - coef(model) %*% x) - link(zeta[z] - coef(model) %*% x))), c(2,3), function(x) weighted.mean(x, wi))
        l3[1:n.x1,2:(n.q+1),1:n.x2,z] <- aperm(apply(l2[,,,z], c(2,3), function(x) quantile(x, probs=quantiles)), c(2,1,3))
        l3[nrow(l3),1,1:n.x2,z] <- l3[n.x1,1,1:n.x2,z] - l3[1,1,1:n.x2,z]
        l3[nrow(l3),2:(n.q+1),1:n.x2,z] <- apply(l2[,n.x1,,z] - l2[,1,,z], 2, function(x) quantile(x, probs=quantiles))
      }
      
      dimnames(l3) <- list(paste(c(rep(paste(x1name," ="),n.x1),paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),c(x1vals,"")),
                           c("est",quantiles),
                           paste(c(rep(paste(x2name,"="),n.x2)),x2vals),
                           paste(c(rep("Y =",n.y)), c(1:n.y))) 
      
      if (is.null(did)){did <- c(x2vals[1],x2vals[n.x2])} else{did <- did}
      l4 <- array(NA, c(n.y,n.q+1))
      for (i in 1:n.y){
        l4[i,1] <- (l3[n.x1,1,match(did[2],x2vals),i] - l3[1,1,match(did[2],x2vals),i]) -  (l3[n.x1,1,match(did[1],x2vals),i] - l3[1,1,match(did[1],x2vals),i])
        l4[i,2:(n.q+1)] <- quantile((l2[ ,n.x1,match(did[2],x2vals),i] - l2[ ,1,match(did[2],x2vals),i]) -  (l2[ ,n.x1,match(did[1],x2vals),i] - l2[ ,1,match(did[1],x2vals),i]), probs=quantiles)
      }
      yvals <- 1:n.y
      dimnames(l4) <- list(paste(c(rep(paste("Y","="),n.y)),yvals),c("est",quantiles)) 
      
      ans <- new("post", 
                 est=round(l3, digits=digits), 
                 did=round(l4, digits=digits), 
                 sims=l2, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
  }
  
  else{
    
    if (is.null(x1name)){
      
      X_temp <- array(NA, c(n.obs,k))
      X <- array(NA, c(n.obs,k-1))
      
      newdata <- data.frame(model$model)
      if (!is.null(holds)){
        for (j in 1:length(holds)){
          newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
        }
      }
      X_temp[ , ] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
      X[ , ] <- X_temp[,-1]
      X <- aperm(X)
      
      l1 <- apply(link(-tau[,cut+1] + beta %*% X), 1, function(x) weighted.mean(x, wi))
      l2 <- array(NA, c(1,n.q+1))
      l2[1,1] <- weighted.mean(link(-zeta[cut+1] + coef(model) %*% X))
      l2[1,2:(n.q+1)] <- quantile(l1, probs=quantiles)
      colnames(l2) <- c("est",quantiles)
      
      ans <- new("post", 
                 est=round(l2, digits=digits), 
                 did=NULL, 
                 sims=l1, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
    
    
    else if (is.null(x2name)){
      
      n.x1 <- length(x1vals)
      
      X_temp <- array(NA, c(n.obs,k,n.x1))
      X <- array(NA, c(n.obs,k-1,n.x1))
      
      for (i in 1:(n.x1)){
        newdata <- data.frame(model$model)
        if (!is.null(holds)){
          for (j in 1:length(holds)){
            newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
          }
        }
        newdata[ ,x1name] <- x1vals[i]
        X_temp[ , ,i] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
        X[ , ,i] <- X_temp[,-1,i]
      }
      
      X <- aperm(X, c(2,1,3))
      l1 <- apply(apply(X, c(2,3), function(x) link(-tau[,cut+1] + beta %*% x)), c(1,3), function(x) weighted.mean(x, wi))
      l2 <- array(NA, c(n.x1+1,n.q+1))
      l2[1:n.x1,1] <- apply(apply(X, c(2,3), function(x) link(-zeta[cut+1] + coef(model) %*% x)), 2, function(x) weighted.mean(x,wi))
      l2[1:n.x1,2:(n.q+1)] <- t(apply(l1, 2, function(x) quantile(x, probs=quantiles)))
      l2[nrow(l2),1] <- l2[n.x1,1] - l2[1,1]
      l2[nrow(l2),2:(n.q+1)] <- quantile(l1[ ,ncol(l1)] - l1[ ,1], probs=quantiles)
      rownames(l2) <- c(paste(c(rep(paste(x1name,"="),n.x1),
                                paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),
                              c(x1vals,"")))  
      colnames(l2) <- c("est",quantiles)
      
      ans <- new("post", 
                 est=round(l2, digits=digits), 
                 did=NULL, 
                 sims=l1, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    } 
    
    else{
      
      n.x1 <- length(x1vals)
      n.x2 <- length(x2vals)
      
      X_temp <- array(NA, c(n.obs,k,n.x1,n.x2))
      X <- array(NA, c(n.obs,k-1,n.x1,n.x2))
      
      for (j in 1:n.x2){
        for (i in 1:n.x1){
          newdata <- data.frame(model$model)
          if (!is.null(holds)){
            for (k in 1:length(holds)){
              newdata[ ,names(holds)[k]] <- as.numeric(holds[k])
            }
          }
          newdata[ ,x1name] <- x1vals[i]
          newdata[ ,x2name] <- x2vals[j]
          X_temp[ , ,i,j] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
          X[ , ,i,j] <- X_temp[,-1,i,j]
        }
      }
      
      X <- aperm(X, c(2,1,3,4))
      l1 <- apply(apply(X, c(2,3,4), function(x) link(-tau[,cut+1] + beta %*% x)), c(1,3,4), function(x) weighted.mean(x, wi))
      l2 <- array(NA, c(n.x1+1,n.q+1,n.x2))
      
      l2[1:n.x1,1,1:n.x2] <- apply(apply(X, c(2,3,4), function(x) link(-zeta[cut+1] + coef(model) %*% x)), c(2,3), function(x) weighted.mean(x, wi))
      l2[1:n.x1,2:(n.q+1),1:n.x2] <- aperm(apply(l1, c(2,3), function(x) quantile(x, probs=quantiles)), c(2,1,3))
      l2[nrow(l2),1,1:n.x2] <- l2[n.x1,1,1:n.x2] - l2[1,1,1:n.x2]
      l2[nrow(l2),2:(n.q+1),1:n.x2] <- apply(l1[ ,n.x1,1:n.x2] - l1[ ,1,1:n.x2], 2, function(x) quantile(x, probs=quantiles))
      
      dimnames(l2) <- list(paste(c(rep(paste(x1name," ="),n.x1),paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),c(x1vals,"")),
                           c("est",quantiles),
                           paste(c(rep(paste(x2name," ="),n.x2)),
                                 c(x2vals)))
      
      if (is.null(did)){did <- c(x2vals[1],x2vals[n.x2])} else{did <- did}
      l3 <- array(NA, c(1,n.q+1))
      l3[1,1] <- (l2[n.x1,1,match(did[2],x2vals)] - l2[1,1,match(did[2],x2vals)]) - (l2[n.x1,1,match(did[1],x2vals)] - l2[1,1,match(did[1],x2vals)])
      l3[1,2:(n.q+1)] <- quantile((l1[ ,n.x1,match(did[2],x2vals)] - l1[ ,1,match(did[2],x2vals)]) -  (l1[ ,n.x1,match(did[1],x2vals)] - l1[ ,1,match(did[1],x2vals)]), probs=quantiles)
      dimnames(l3) <- list("did",c("est",quantiles)) 
      
      ans <- new("post", 
                 est=round(l2, digits=digits), 
                 did=round(l3, digits=digits), 
                 sims=l1, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
  }
}

setMethod("post", signature(model = "lm"), post.glm)
setMethod("post", signature(model = "glm"), post.glm)
setMethod("post", signature(model = "svyglm"), post.glm)  
setMethod("post", signature(model = "polr"), post.polr)  

#### postSim FUNCTION ####

setGeneric("postSim", 
           function(object, n.sims=1000){
             standardGeneric("postSim")
           }
)

setClass("postSim", 
         slots = c(coef = "matrix",
                   sigma = "numeric")
)

setClass("postSim.polr", 
         slots = c(coef = "matrix",
                   zeta = "matrix")
)

setMethod("postSim", signature(object = "lm"),
          function(object, n.sims=1000)
          {
            object.class <- class(object)[[1]]
            summ <- summary (object)
            coef <- summ$coef[,1:2,drop=FALSE]
            dimnames(coef)[[2]] <- c("coef.est","coef.sd")
            sigma.hat <- summ$sigma
            beta.hat <- coef[,1,drop = FALSE]
            V.beta <- summ$cov.unscaled
            n <- summ$df[1] + summ$df[2]
            k <- summ$df[1]
            sigma <- rep (NA, n.sims)
            beta <- array (NA, c(n.sims,k))
            dimnames(beta) <- list (NULL, rownames(beta.hat))
            for (s in 1:n.sims){
              sigma[s] <- sigma.hat*sqrt((n-k)/rchisq(1,n-k))
              beta[s,] <- MASS::mvrnorm (1, beta.hat, V.beta*sigma[s]^2)
            }
            ans <- new("postSim",
                       coef = beta,
                       sigma = sigma)
            return (ans)
          }
)

setMethod("postSim", signature(object = "glm"),
          function(object, n.sims=1000)
          {
            object.class <- class(object)[[1]]
            summ <- summary (object, correlation=TRUE, dispersion = object$dispersion)
            coef <- summ$coef[,1:2,drop=FALSE]
            dimnames(coef)[[2]] <- c("coef.est","coef.sd")
            beta.hat <- coef[,1,drop=FALSE]
            sd.beta <- coef[,2,drop=FALSE]
            corr.beta <- summ$corr
            n <- summ$df[1] + summ$df[2]
            k <- summ$df[1]
            V.beta <- corr.beta * array(sd.beta,c(k,k)) * t(array(sd.beta,c(k,k)))
            beta <- MASS::mvrnorm (n.sims, beta.hat, V.beta)
            beta2 <- array (0, c(n.sims,length(coefficients(object))))
            dimnames(beta2) <- list (NULL, names(coefficients(object)))
            beta2[,dimnames(beta2)[[2]]%in%dimnames(beta)[[2]]] <- beta
            sigma <- rep (sqrt(summ$dispersion), n.sims)
            ans <- new("postSim",
                       coef = beta2,
                       sigma = sigma)
            return(ans)
          }
)

setMethod("postSim", signature(object = "polr"),
          function(object, n.sims=1000){
            x <- as.matrix(model.matrix(object))
            coefs <- coef(object)
            k <- length(coefs)
            zeta <- object$zeta
            Sigma <- vcov(object)
            
            if(n.sims==1){
              parameters <- t(MASS::mvrnorm(n.sims, c(coefs, zeta), Sigma))
            }else{
              parameters <- MASS::mvrnorm(n.sims, c(coefs, zeta), Sigma)
            }
            ans <- new("postSim.polr",
                       coef = parameters[,1:k,drop=FALSE],
                       zeta = parameters[,-(1:k),drop=FALSE])
            return(ans)
          }
)


setMethod("postSim", signature(object = "svyglm"),
          function(object, n.sims=1000)
          {
            object.class <- class(object)[[1]]
            if (family(object)[1]=="gaussian"){
              summ <- summary (object)
              coef <- summ$coef[,1:2,drop=FALSE]
              dimnames(coef)[[2]] <- c("coef.est","coef.sd")
              sigma.hat <- sigma(object)
              beta.hat <- coef[,1,drop = FALSE]
              V.beta <- (summ$cov.scaled)*(1/(sigma.hat)^2) #$cov.unscaled is actually scaled in svyglm!
              n <- summ$df[1] + summ$df[2]
              k <- summ$df[1]
              sigma <- rep (NA, n.sims)
              beta <- array (NA, c(n.sims,k))
              dimnames(beta) <- list (NULL, rownames(beta.hat))
              for (s in 1:n.sims){
                sigma[s] <- sigma.hat*sqrt((n-k)/rchisq(1,n-k))
                beta[s,] <- MASS::mvrnorm (1, beta.hat, V.beta*sigma[s]^2)
              }
              ans <- new("postSim",
                         coef = beta,
                         sigma = sigma)
              return (ans)
            }
            else{
              summ <- summary (object, correlation=TRUE, dispersion = object$dispersion)
              coef <- summ$coef[,1:2,drop=FALSE]
              dimnames(coef)[[2]] <- c("coef.est","coef.sd")
              beta.hat <- coef[,1,drop=FALSE]
              sd.beta <- coef[,2,drop=FALSE]
              corr.beta <- summ$corr
              n <- summ$df[1] + summ$df[2]
              k <- summ$df[1]
              V.beta <- corr.beta * array(sd.beta,c(k,k)) * t(array(sd.beta,c(k,k)))
              beta <- MASS::mvrnorm (n.sims, beta.hat, V.beta)
              beta2 <- array (0, c(n.sims,length(coefficients(object))))
              dimnames(beta2) <- list (NULL, names(coefficients(object)))
              beta2[,dimnames(beta2)[[2]]%in%dimnames(beta)[[2]]] <- beta
              sigma <- rep (sqrt(summ$dispersion), n.sims)
              
              ans <- new("postSim",
                         coef = beta2,
                         sigma = sigma)
              return(ans)
            }
          }
)


#### ANES-GSS 2020 ANALYSIS ####

## Create Survey Design to Weight Data
ANESGSS_Design <- svydesign(id=~1, 
                         weights=~weight, 
                         nest=TRUE, 
                         data=ANESGSS2020)

## Regression Models
mod_concern <- svyglm(covidconcern ~ auth_scale + age + male + black + hispanic + education + income + unemp + south + engagement + auth_scale*engagement, design = ANESGSS_Design)
mod2_leftright_ANESGSS <- svyglm(dv_leftright ~ auth_scale + age + male + black + hispanic + education +  income + unemp + south + engagement + auth_scale*engagement, design = ANESGSS_Design)

## Conditional Marginal Effects
concern_ANESGSS <- post(mod_concern, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)
leftright_ANESGSS <- post(mod2_leftright_ANESGSS, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)


## Concern over COVID-19 (FIGURE 2)
tiff(file="concern_ANESGSS.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(concern_ANESGSS@est[3,1,1], concern_ANESGSS@est[3,1,2], concern_ANESGSS@est[3,1,3], concern_ANESGSS@est[3,1,4], concern_ANESGSS@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.49,.45), 
     xlim=c(0,1))
axis(2, at=c(-.6,-.4,-.2,0,.2,.4, .6), labels=c("","-.40","-.20",".00",".20",".40", ""), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(concern_ANESGSS@est[3,2,1], concern_ANESGSS@est[3,2,2], concern_ANESGSS@est[3,2,3], concern_ANESGSS@est[3,2,4], concern_ANESGSS@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(concern_ANESGSS@est[3,3,1], concern_ANESGSS@est[3,3,2], concern_ANESGSS@est[3,3,3], concern_ANESGSS@est[3,3,4], concern_ANESGSS@est[3,3,5]))
abline(h=0, lty=3)
title("ANES-GSS 2020", line = -.1)

dev.off()

## Left-Right Orientation (FIGURE 1)
tiff(file="leftright_ANESGSS.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(leftright_ANESGSS@est[3,1,1], leftright_ANESGSS@est[3,1,2], leftright_ANESGSS@est[3,1,3], leftright_ANESGSS@est[3,1,4], leftright_ANESGSS@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.60), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4, .6), labels=c("-.40","-.20",".00",".20",".40", ".60"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(leftright_ANESGSS@est[3,2,1], leftright_ANESGSS@est[3,2,2], leftright_ANESGSS@est[3,2,3], leftright_ANESGSS@est[3,2,4], leftright_ANESGSS@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(leftright_ANESGSS@est[3,3,1], leftright_ANESGSS@est[3,3,2], leftright_ANESGSS@est[3,3,3], leftright_ANESGSS@est[3,3,4], leftright_ANESGSS@est[3,3,5]))
abline(h=0, lty=3)
title("ANES-GSS 2020", line = -.1)

dev.off()

#### LUCID 2020 ANALYSIS ####

## Regression Models
mod2_behavior_2020 <- lm(behavior_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2020)
mod2_econ_2020 <- lm(econ_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2020)
mod2_publichealth_2020 <- lm(publichealth_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2020)
mod2_leftright_2020 <- lm(dv_leftright ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2020)

## Conditional Marginal Effects
behavior_2020 <- post(mod2_behavior_2020, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)
econ_2020 <- post(mod2_econ_2020, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)
publichealth_2020 <- post(mod2_publichealth_2020, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)
leftright_2020 <- post(mod2_leftright_2020, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)

## Health Behaviors (FIGURE 3)
tiff(file="behavior_2020.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(behavior_2020@est[3,1,1], behavior_2020@est[3,1,2], behavior_2020@est[3,1,3], behavior_2020@est[3,1,4], behavior_2020@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.4), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4), labels=c("-.40","-.20",".00",".20",".40"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(behavior_2020@est[3,2,1], behavior_2020@est[3,2,2], behavior_2020@est[3,2,3], behavior_2020@est[3,2,4], behavior_2020@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(behavior_2020@est[3,3,1], behavior_2020@est[3,3,2], behavior_2020@est[3,3,3], behavior_2020@est[3,3,4], behavior_2020@est[3,3,5]))
abline(h=0, lty=3)
title("Masking and Distancing (2020)", line = -.1)

dev.off()

## Economic Interventions (FIGURE 5)
tiff(file="econ_2020.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(econ_2020@est[3,1,1], econ_2020@est[3,1,2], econ_2020@est[3,1,3], econ_2020@est[3,1,4], econ_2020@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.4), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4), labels=c("-.40","-.20",".00",".20",".40"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(econ_2020@est[3,2,1], econ_2020@est[3,2,2], econ_2020@est[3,2,3], econ_2020@est[3,2,4], econ_2020@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(econ_2020@est[3,3,1], econ_2020@est[3,3,2], econ_2020@est[3,3,3], econ_2020@est[3,3,4], econ_2020@est[3,3,5]))
abline(h=0, lty=3)
title("Lucid 2020", line = -.1)

dev.off()

## Public Health Measures (FIGURE 4)
tiff(file="publichealth_2020.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(publichealth_2020@est[3,1,1], publichealth_2020@est[3,1,2], publichealth_2020@est[3,1,3], publichealth_2020@est[3,1,4], publichealth_2020@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.4), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4), labels=c("-.40","-.20",".00",".20",".40"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(publichealth_2020@est[3,2,1], publichealth_2020@est[3,2,2], publichealth_2020@est[3,2,3], publichealth_2020@est[3,2,4], publichealth_2020@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(publichealth_2020@est[3,3,1], publichealth_2020@est[3,3,2], publichealth_2020@est[3,3,3], publichealth_2020@est[3,3,4], publichealth_2020@est[3,3,5]))
abline(h=0, lty=3)
title("Lucid 2020", line = -.1)

dev.off()

## Left-Right Orientation (FIGURE 1)
tiff(file="leftright_2020.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(leftright_2020@est[3,1,1], leftright_2020@est[3,1,2], leftright_2020@est[3,1,3], leftright_2020@est[3,1,4], leftright_2020@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.60), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4,.6), labels=c("-.40","-.20",".00",".20",".40", ".60"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(leftright_2020@est[3,2,1], leftright_2020@est[3,2,2], leftright_2020@est[3,2,3], leftright_2020@est[3,2,4], leftright_2020@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(leftright_2020@est[3,3,1], leftright_2020@est[3,3,2], leftright_2020@est[3,3,3], leftright_2020@est[3,3,4], leftright_2020@est[3,3,5]))
abline(h=0, lty=3)
title("Lucid 2020", line = -.1)

dev.off()


#### LUCID 2021 ANALYSIS ####

## Create New Dataset without missing values for Logistic Model
Lucid2021_Vaccination <- na.omit(subset(Lucid2021[, c("dv_vaccinated", "auth_scale", "age01",
                                                      "male", "black", "hispanic", "educ01",
                                                      "income01", "unemp", "south", "covid_posandsus", 
                                                      "engagement")]))

## Regression Models
mod2_behavior_2021 <- lm(behavior_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement, data = Lucid2021)
mod2_vaccination_2021 <- glm(as.factor(dv_vaccinated) ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement, data = Lucid2021_Vaccination, family = binomial(link = "logit"))
mod2_econ_2021 <- lm(econ_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement, data = Lucid2021)
mod2_publichealth_2021 <- lm(publichealth_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement, data = Lucid2021)
mod2_leftright_2021 <- lm(dv_leftright ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2021)

## Conditional Marginal Effects
behavior_2021 <- post(mod2_behavior_2021, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)
econ_2021 <- post(mod2_econ_2021, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)
publichealth_2021 <- post(mod2_publichealth_2021, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)
leftright_2021 <- post(mod2_leftright_2021, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)
vaccination_2021 <- post(mod2_vaccination_2021, x1name =  "auth_scale", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)

## Health Behaviors (FIGURE 3)
tiff(file="behavior_2021.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(behavior_2021@est[3,1,1], behavior_2021@est[3,1,2], behavior_2021@est[3,1,3], behavior_2021@est[3,1,4], behavior_2021@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.4), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4), labels=c("-.40","-.20",".00",".20",".40"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(behavior_2021@est[3,2,1], behavior_2021@est[3,2,2], behavior_2021@est[3,2,3], behavior_2021@est[3,2,4], behavior_2021@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(behavior_2021@est[3,3,1], behavior_2021@est[3,3,2], behavior_2021@est[3,3,3], behavior_2021@est[3,3,4], behavior_2021@est[3,3,5]))
abline(h=0, lty=3)
title("Masking and Distancing (2021)", line = -.1)

dev.off()

## Vaccination (FIGURE 3)
tiff(file="vaccination_2021.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(vaccination_2021@est[3,1,1], vaccination_2021@est[3,1,2], vaccination_2021@est[3,1,3], vaccination_2021@est[3,1,4], vaccination_2021@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.4), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4), labels=c("-.40","-.20",".00",".20",".40"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(vaccination_2021@est[3,2,1], vaccination_2021@est[3,2,2], vaccination_2021@est[3,2,3], vaccination_2021@est[3,2,4], vaccination_2021@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(vaccination_2021@est[3,3,1], vaccination_2021@est[3,3,2], vaccination_2021@est[3,3,3], vaccination_2021@est[3,3,4], vaccination_2021@est[3,3,5]))
abline(h=0, lty=3)
title("COVID-19 Vaccination (2021)", line = -.1)

dev.off()

## Economic Interventions (FIGURE 5)
tiff(file="econ_2021.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(econ_2021@est[3,1,1], econ_2021@est[3,1,2], econ_2021@est[3,1,3], econ_2021@est[3,1,4], econ_2021@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.4), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4), labels=c("-.40","-.20",".00",".20",".40"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(econ_2021@est[3,2,1], econ_2021@est[3,2,2], econ_2021@est[3,2,3], econ_2021@est[3,2,4], econ_2021@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(econ_2021@est[3,3,1], econ_2021@est[3,3,2], econ_2021@est[3,3,3], econ_2021@est[3,3,4], econ_2021@est[3,3,5]))
abline(h=0, lty=3)
title("Lucid 2021", line = -.1)

dev.off()

## Public Health Measures (FIGURE 4)
tiff(file="publichealth_2021.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(publichealth_2021@est[3,1,1], publichealth_2021@est[3,1,2], publichealth_2021@est[3,1,3], publichealth_2021@est[3,1,4], publichealth_2021@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.4), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4), labels=c("-.40","-.20",".00",".20",".40"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(publichealth_2021@est[3,2,1], publichealth_2021@est[3,2,2], publichealth_2021@est[3,2,3], publichealth_2021@est[3,2,4], publichealth_2021@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(publichealth_2021@est[3,3,1], publichealth_2021@est[3,3,2], publichealth_2021@est[3,3,3], publichealth_2021@est[3,3,4], publichealth_2021@est[3,3,5]))
abline(h=0, lty=3)
title("Lucid 2021", line = -.1)

dev.off()

## Left-Right Orientation (FIGURE 1)
tiff(file="leftright_2021.tiff", type="cairo", height = 5, width = 5, family="serif", units = 'in', res=1500, pointsize=11)

plot(x = c(0,.25,.5,.75,1), 
     y=c(leftright_2021@est[3,1,1], leftright_2021@est[3,1,2], leftright_2021@est[3,1,3], leftright_2021@est[3,1,4], leftright_2021@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Engagement", 
     axes=FALSE, 
     ylim=c(-.4,.60), 
     xlim=c(0,1))
axis(2, at=c(-.4,-.2,0,.2,.4,.6), labels=c("-.40","-.20",".00",".20",".40",".60"), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(leftright_2021@est[3,2,1], leftright_2021@est[3,2,2], leftright_2021@est[3,2,3], leftright_2021@est[3,2,4], leftright_2021@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(leftright_2021@est[3,3,1], leftright_2021@est[3,3,2], leftright_2021@est[3,3,3], leftright_2021@est[3,3,4], leftright_2021@est[3,3,5]))
abline(h=0, lty=3)
title("Lucid 2021", line = -.1)

dev.off()


#### APPENDIX 3: REGRESSION TABLES ####

# Rename ANES-GSS covariates so they are identical across datasets
ANESGSS2020$age01 <- ANESGSS2020$age
ANESGSS2020$educ01 <- ANESGSS2020$education
ANESGSS2020$income01 <- ANESGSS2020$income
ANESGSS_Design <- svydesign(id=~1, 
                            weights=~weight, 
                            nest=TRUE, 
                            data=ANESGSS2020)

## Regressions
mod_concern_2020 <- svyglm(covidconcern ~ auth_scale + age01 + male + black + hispanic + educ01 + income01 + unemp + south + engagement + auth_scale*engagement, design = ANESGSS_Design)
mod2_behavior_2020 <- lm(behavior_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2020)
mod2_econ_2020 <- lm(econ_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2020)
mod2_publichealth_2020 <- lm(publichealth_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2020)
mod2_behavior_2021 <- lm(behavior_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement, data = Lucid2021)
mod2_vaccination_2021 <- glm(as.factor(dv_vaccinated) ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement, data = Lucid2021, family = binomial(link = "logit"))
mod2_econ_2021 <- lm(econ_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement, data = Lucid2021)
mod2_publichealth_2021 <- lm(publichealth_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement, data = Lucid2021)
mod2_leftright_ANESGSS <- svyglm(dv_leftright ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, design = ANESGSS_Design)
mod2_leftright_2020 <- lm(dv_leftright ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2020)
mod2_leftright_2021 <- lm(dv_leftright ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement, data = Lucid2021)

## Regression Tables 
mod2_2020 <- stargazer(mod_concern_2020, mod2_behavior_2020, mod2_publichealth_2020, mod2_econ_2020,
                       title="",
                       dep.var.labels=c("COVID Concern","Health Behaviors", "Public Health Measures", "Economic Interventions"),
                       covariate.labels=c("Authoritarianism", "Age", "Male", "Black", "Hispanic", 
                                          "Education", "Income", "Unemployed", "South", "Engagement",
                                          "Authoritarianism:Engagement"),
                       keep.stat=c("n"),
                       report = "vcsp", 
                       digits = 3,
                       digits.extra = 0,
                       align=TRUE, 
                       no.space=TRUE,
                       type = 'html',
                       header = FALSE, 
                       out = 'table_2020.html')

mod2_2021 <- stargazer(mod2_behavior_2021, mod2_vaccination_2021, mod2_publichealth_2021, mod2_econ_2021,
                       title="",
                       dep.var.labels=c("Health Behaviors", "Vaccination Status", "Public Health Measures", "Economic Interventions"),
                       covariate.labels=c("Authoritarianism", "Age", "Male", "Black", "Hispanic", 
                                          "Education", "Income", "Unemployed", "South", "Prior COVID Case",
                                          "Engagement", "Authoritarianism:Engagement"),
                       keep.stat=c("n"),
                       report = "vcsp", 
                       digits = 3,
                       digits.extra = 0,
                       align=TRUE, 
                       no.space=TRUE,
                       type = 'html',
                       header = FALSE, 
                       out = 'table_2021.html')

mod2_lr <- stargazer(mod2_leftright_ANESGSS, mod2_leftright_2020, mod2_leftright_2021,
                     title="",
                     dep.var.labels=c("ANES-GSS","Lucid 2020", "Lucid 2021"),
                     covariate.labels=c("Authoritarianism", "Age", "Male", "Black", "Hispanic", 
                                        "Education", "Income", "Unemployed", "South", "Engagement",
                                        "Authoritarianism:Engagement"),
                     keep.stat=c("n"),
                     report = "vcsp", 
                     digits = 3,
                     digits.extra = 0,
                     align=TRUE, 
                     no.space=TRUE,
                     type = 'html',
                     header = FALSE, 
                     out = 'table_lr.html')


#### APPENDIX 4: ANALYSIS WITH PID/IDEO CONTROLS ####

## Consistent covariate names
ANESGSS2020$engagement <- ANESGSS2020$knowledge
ANESGSS2020$age01 <- ANESGSS2020$age
ANESGSS2020$educ01 <- ANESGSS2020$education
ANESGSS2020$income01 <- ANESGSS2020$income
ANESGSS2020$pid01 <- ANESGSS2020$pid01_t2
ANESGSS2020$ideo <- ANESGSS2020$ideo01_t2
ANESGSS_Design <- svydesign(id=~1, 
                            weights=~weight, 
                            nest=TRUE, 
                            data=ANESGSS2020)

## Regression Models
mod_concern_lr <- svyglm(covidconcern ~ auth_scale + age01 + male + black + hispanic + educ01 + income01 + unemp + south + engagement + auth_scale*engagement + pid01 + ideo, design = ANESGSS_Design)
mod2_behavior_2020_lr <- lm(behavior_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement + pid01 + ideo, data = Lucid2020)
mod2_econ_2020_lr <- lm(econ_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement + pid01 + ideo, data = Lucid2020)
mod2_publichealth_2020_lr <- lm(publichealth_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement + auth_scale*engagement + pid01 + ideo, data = Lucid2020)
mod2_behavior_2021_lr <- lm(behavior_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement + pid01 + ideo, data = Lucid2021)
mod2_vaccination_2021_lr <- glm(as.factor(dv_vaccinated) ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement + pid01 + ideo, data = Lucid2021, family = binomial(link = "logit"))
mod2_econ_2021_lr <- lm(econ_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement + pid01 + ideo, data = Lucid2021)
mod2_publichealth_2021_lr <- lm(publichealth_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement + auth_scale*engagement + pid01 + ideo, data = Lucid2021)


##### APPENDIX 5: BINNING ESTIMATEs #####
ANESGSS2020 <- as.data.frame(ANESGSS2020)

## Left-Right Orientation
lr_anesgss <- interflex(data = ANESGSS2020, Y = "dv_leftright", D = "auth_scale", X = "engagement", weights = "weight",
                        Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south"), 
                        estimator = "linear",  vcov.type = "robust", na.rm = T,
                        theme.bw = TRUE, show.grid = FALSE, 
                        main = "", 
                        xlab = "Political Engagement", ylab = "Marginal Effect")

plot(lr_anesgss$figure)
dev.off()

## Concern over COVID-19
concern_2020 <- interflex(data = ANESGSS2020, Y = "covidconcern", D = "auth_scale", X = "engagement", weights = "weight",
                          Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south"), 
                          estimator = "linear",  vcov.type = "robust", na.rm = T,
                          theme.bw = TRUE, show.grid = FALSE, 
                          main = "", 
                          xlab = "Political Engagement", ylab = "Marginal Effect")

plot(concern_2020$figure)
dev.off()

## Health Behavior
behavior_2020 <- interflex(data = Lucid2020, Y = "behavior_scale", D = "auth_scale", X = "engagement", 
                           Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south"), 
                           estimator = "binning", nbins = 3, vcov.type = "robust", na.rm = T,
                           theme.bw = TRUE, show.grid = FALSE, 
                           main = "Masking and Distancing (2020)", 
                           xlab = "Political Engagement", ylab = "Marginal Effect")

plot(behavior_2020$figure)
dev.off()

## Public Health
publichealth_2020 <- interflex(data = Lucid2020, Y = "publichealth_scale", D = "auth_scale", X = "engagement", 
                               Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south"), 
                               estimator = "binning", nbins = 3, vcov.type = "robust", na.rm = T,
                               theme.bw = TRUE, show.grid = FALSE, 
                               main = "Lucid 2020", 
                               xlab = "Political Engagement", ylab = "Marginal Effect")

plot(publichealth_2020$figure)
dev.off()

## Economic Interventions
econ_2020 <- interflex(data = Lucid2020, Y = "econ_scale", D = "auth_scale", X = "engagement", 
                       Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south"), 
                       estimator = "binning", nbins = 3, vcov.type = "robust", na.rm = T,
                       theme.bw = TRUE, show.grid = FALSE, 
                       main = "Lucid 2020", 
                       xlab = "Political Engagement", ylab = "Marginal Effect")

plot(econ_2020$figure)
dev.off()

## Left-Right Orientation
lr_2020 <- interflex(data = Lucid2020, Y = "dv_leftright", D = "auth_scale", X = "engagement", 
                     Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south"), 
                     estimator = "binning", nbins = 3, vcov.type = "robust", na.rm = T,
                     theme.bw = TRUE, show.grid = FALSE, 
                     main = "Lucid 2020", 
                     xlab = "Political Engagement", ylab = "Marginal Effect")

plot(lr_2020$figure)
dev.off()

## Health Behavior
behavior_2021 <- interflex(data = Lucid2021, Y = "behavior_scale", D = "auth_scale", X = "engagement", 
                           Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south", "covid_posandsus"), 
                           estimator = "binning", nbins = 3, vcov.type = "robust", na.rm = T,
                           theme.bw = TRUE, show.grid = FALSE, 
                           main = "Masking and Distancing (2021)", 
                           xlab = "Political Engagement", ylab = "Marginal Effect")

plot(behavior_2021$figure)
dev.off()

## Public Health
publichealth_2021 <- interflex(data = Lucid2021, Y = "publichealth_scale", D = "auth_scale", X = "engagement", 
                               Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south", "covid_posandsus"), 
                               estimator = "binning", nbins = 3, vcov.type = "robust", na.rm = T,
                               theme.bw = TRUE, show.grid = FALSE, 
                               main = "Lucid 2021", 
                               xlab = "Political Engagement", ylab = "Marginal Effect")

plot(publichealth_2021$figure)
dev.off()

## Economic Interventions
econ_2021 <- interflex(data = Lucid2021, Y = "econ_scale", D = "auth_scale", X = "engagement", 
                       Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south", "covid_posandsus"), 
                       estimator = "binning", nbins = 3, vcov.type = "robust", na.rm = T,
                       theme.bw = TRUE, show.grid = FALSE, 
                       main = "Lucid 2021", 
                       xlab = "Political Engagement", ylab = "Marginal Effect")

plot(econ_2021$figure)
dev.off()

## Vaccination Status
vaccination_2021 <- interflex(data = Lucid2021, Y = "dv_vaccinated", D = "auth_scale", X = "engagement", 
                              Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south", "covid_posandsus"), 
                              estimator = "binning", method="logit", nbins = 3, vcov.type = "robust", na.rm = T,
                              theme.bw = TRUE, show.grid = FALSE, 
                              main = "COVID-19 Vaccination (2021)", 
                              xlab = "Political Engagement", ylab = "Marginal Effect")

plot(vaccination_2021$figure)
dev.off()

## Left-Right Orientation
lr_2021 <- interflex(data = Lucid2021, Y = "dv_leftright", D = "auth_scale", X = "engagement", 
                     Z = c("age01", "male", "black", "hispanic", "educ01", "income01", "unemp", "south"), 
                     estimator = "binning", nbins = 3, vcov.type = "robust", na.rm = T,
                     theme.bw = TRUE, show.grid = FALSE, 
                     main = "Lucid 2021", 
                     xlab = "Political Engagement", ylab = "Marginal Effect")

plot(lr_2021$figure)
dev.off()


#### APPENDIX 8: NO-INTERACTION ANALYSIS ####

## Regression Models
mod1_concern_2020 <- svyglm(covidconcern ~ auth_scale + age + male + black + hispanic + education + income + unemp + south + engagement, design = ANESGSS_Design)
mod1_behavior_2020 <- lm(behavior_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement, data = Lucid2020)
mod1_econ_2020 <- lm(econ_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement, data = Lucid2020)
mod1_publichealth_2020 <- lm(publichealth_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + engagement, data = Lucid2020)
mod1_behavior_2021 <- lm(behavior_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement, data = Lucid2021)
mod1_vaccination_2021 <- lm(dv_vaccinated ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement, data = Lucid2021)
mod1_econ_2021 <- lm(econ_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement, data = Lucid2021)
mod1_publichealth_2021 <- lm(publichealth_scale ~ auth_scale + age01 + male + black + hispanic + educ01 +  income01 + unemp + south + covid_posandsus + engagement, data = Lucid2021)

## Put Parameters/Errors in Data Frame
model1Frame <- data.frame(Variable = rownames(summary(mod1_behavior_2020)$coef),
                          Coefficient = summary(mod1_behavior_2020)$coef[, 1],
                          SE = summary(mod1_behavior_2020)$coef[, 2],
                          Domain = "Health Behaviors")
model1Frame <- model1Frame[2,]

model2Frame <- data.frame(Variable = rownames(summary(mod1_publichealth_2020)$coef),
                          Coefficient = summary(mod1_publichealth_2020)$coef[, 1],
                          SE = summary(mod1_publichealth_2020)$coef[, 2],
                          Domain = "Public Health Measures")
model2Frame <- model2Frame[2,]

model3Frame <- data.frame(Variable = rownames(summary(mod1_econ_2020)$coef),
                          Coefficient = summary(mod1_econ_2020)$coef[, 1],
                          SE = summary(mod1_econ_2020)$coef[, 2],
                          Domain = "Economic Interventions")
model3Frame <- model3Frame[2,]

model4Frame <- data.frame(Variable = rownames(summary(mod1_concern_2020)$coef),
                          Coefficient = summary(mod1_concern_2020)$coef[, 1],
                          SE = summary(mod1_concern_2020)$coef[, 2],
                          Domain = "COVID-19 Concern")
model4Frame <- model4Frame[2,]

model5Frame <- data.frame(Variable = rownames(summary(mod1_behavior_2021)$coef),
                          Coefficient = summary(mod1_behavior_2021)$coef[, 1],
                          SE = summary(mod1_behavior_2021)$coef[, 2],
                          Domain = "Health Behaviors")
model5Frame <- model5Frame[2,]

model6Frame <- data.frame(Variable = rownames(summary(mod1_vaccination_2021)$coef),
                          Coefficient = summary(mod1_vaccination_2021)$coef[, 1],
                          SE = summary(mod1_vaccination_2021)$coef[, 2],
                          Domain = "Vaccination Status")
model6Frame <- model6Frame[2,]

model7Frame <- data.frame(Variable = rownames(summary(mod1_publichealth_2021)$coef),
                          Coefficient = summary(mod1_publichealth_2021)$coef[, 1],
                          SE = summary(mod1_publichealth_2021)$coef[, 2],
                          Domain = "Public Health Measures")
model7Frame <- model7Frame[2,]

model8Frame <- data.frame(Variable = rownames(summary(mod1_econ_2021)$coef),
                          Coefficient = summary(mod1_econ_2021)$coef[, 1],
                          SE = summary(mod1_econ_2021)$coef[, 2],
                          Domain = "Economic Interventions")
model8Frame <- model8Frame[2,]

# Figure 8A (2020)
allModelFrame2020 <- data.frame(rbind(model1Frame, model2Frame, model3Frame, model4Frame))
interval2 <- -qnorm((1-0.95)/2) # set 95% CI
zp1 <- ggplot(allModelFrame2020, aes(x = Domain, y = Coefficient), shape = 1)
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp1 <- zp1 + geom_pointrange(aes(x = Domain, y = Coefficient, ymin = Coefficient - SE*interval2, 
                                 ymax = Coefficient + SE*interval2), shape = 1, lwd = .75, position = position_dodge(width = 1/2))
zp1 <- zp1 + theme_bw() + theme(text=element_text(family="serif", face="bold", size=16)) +
  theme(axis.title.x = element_blank())
zp1 <- zp1 + scale_y_continuous(limits = c(-.25,.10), breaks = seq(-.20, .10, by = .10))
zp1 <- zp1 + ggtitle("2020") + theme(plot.title = element_text(hjust = 0.5))
zp1

# Figure 8A (2021)
allModelFrame2021 <- data.frame(rbind(model5Frame, model6Frame, model7Frame, model8Frame))
interval2 <- -qnorm((1-0.95)/2) # set 95% CI
zp2 <- ggplot(allModelFrame2021, aes(x = Domain, y = Coefficient), shape = 1)
zp2 <- zp2 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp2 <- zp2 + geom_pointrange(aes(x = Domain, y = Coefficient, ymin = Coefficient - SE*interval2, 
                                 ymax = Coefficient + SE*interval2), shape = 1, lwd = .75, position = position_dodge(width = 1/2))
zp2 <- zp2 + theme_bw() + theme(text=element_text(family="serif", face="bold", size=16)) +
  theme(axis.title.x = element_blank())
zp2 <- zp2 + scale_y_continuous(limits = c(-.25,.10), breaks = seq(-.20, .20, by = .10)) 
zp2 <- zp2 + ggtitle("2021") + theme(plot.title = element_text(hjust = 0.5))
zp2


